Processing Pipeline for B-cell Data from Popescu et al (2019)

Bacher Group

Author
Affiliation

Jack R. Leary

University of Florida

Published

October 6, 2024

1 Libraries

Before we do anything else we load in the necessary software packages in both R & Python.

1.1 R

Code 
library(dplyr)       # data manipulation
library(Seurat)      # scRNAseq tools
library(scLANE)      # trajectory DE
library(Lamian)      # more trajectory DE
library(ggplot2)     # pretty plots
library(biomaRt)     # gene annotation
library(tradeSeq)    # even more trajectory DE
library(patchwork)   # plot alignment
library(slingshot)   # pseudotime estimation
library(reticulate)  # Python interface
select <- dplyr::select
rename <- dplyr::rename

1.2 Python

Code 
import scvi                                                     # integration
import warnings                                                 # filter out warnings
import numpy as np                                              # linear algebra tools
import scanpy as sc                                             # scRNA-seq tools
import pandas as pd                                             # DataFrames
import anndata as ad                                            # scRNA-seq data structures
import scvelo as scv                                            # RNA velocity
import cellrank as cr                                           # cell fate estimation
import matplotlib.pyplot as plt                                 # pretty plots
from scipy.io import mmread, mmwrite                            # sparse matrix IO
from cellrank.estimators import GPCCA                           # CellRank estimator
from scipy.sparse import coo_matrix, csr_matrix                 # sparse matrices
from cellrank.kernels import PseudotimeKernel, CytoTRACEKernel  # CellRank kernels

Here we set a few global variables to reduce the amount of printed output caused by our Python code.

Code 
warnings.simplefilter('ignore', category=UserWarning)
sc.settings.verbosity = 0
scv.settings.verbosity = 0
cr.settings.verbosity = 0

2 Visualization tools

To make our visualizations prettier we’ll define some consistent color palettes, and set a theme for matplotlib that matches ggplot2::theme_classic().

2.1 Color palettes

Code 
palette_heatmap <- paletteer::paletteer_d("MetBrewer::Cassatt1", direction = -1)
palette_celltype <- as.character(paletteer::paletteer_d("ggsci::category10_d3"))[-c(1:6)]
palette_cluster <- paletteer::paletteer_d("ggsci::default_igv")

2.2 Theme for matplotlib

Code 
base_size = 12
plt.rcParams.update({
    # font
    'font.size': base_size, 
    'font.weight': 'normal',
    # figure
    'figure.dpi': 320, 
    'figure.edgecolor': 'white', 
    'figure.facecolor': 'white', 
    'figure.figsize': (6, 4), 
    'figure.constrained_layout.use': True,
    # axes
    'axes.edgecolor': 'black',
    'axes.grid': False,
    'axes.labelpad': 2.75,
    'axes.labelsize': base_size * 0.8,
    'axes.linewidth': 1.5,
    'axes.spines.right': False,
    'axes.spines.top': False,
    'axes.titlelocation': 'left',
    'axes.titlepad': 11,
    'axes.titlesize': base_size,
    'axes.titleweight': 'normal',
    'axes.xmargin': 0.1, 
    'axes.ymargin': 0.1, 
    # legend
    'legend.borderaxespad': 1,
    'legend.borderpad': 0.5,
    'legend.columnspacing': 2,
    'legend.fontsize': base_size * 0.8,
    'legend.frameon': False,
    'legend.handleheight': 1,
    'legend.handlelength': 1.2,
    'legend.labelspacing': 1,
    'legend.title_fontsize': base_size, 
    'legend.markerscale': 1.5
})

3 Helper functions

We first define a utility function to make our plot legends cleaner to read & easier to make.

Code 
guide_umap <- function(key.size = 4) {
  ggplot2::guides(color = ggplot2::guide_legend(override.aes = list(size = key.size,
                                                                    alpha = 1, 
                                                                    stroke = 0.25)))
}

4 Data

4.1 Import hematopoesis data

We start by importing the AnnData object we downloaded from the Human Developmental Cell Atlas portal; it contains a variety of celltypes beyond just the myeloid lineage which we’re interested.

Code 
ad_liver = ad.read_h5ad('../../data/hematopoesis/fetal_liver_alladata.h5ad')

We subset to just the B-cell lineage, and remove cells with a percentage mitochondrial reads greater than 15%.

Code 
ad_bcell = ad_liver[ad_liver.obs['cell.labels'].isin(['HSC_MPP', 'pre-B cell', 'pro-B cell', 'Pre pro B cell', 'B cell'])]
ad_bcell.obs['celltype'] = (
    ad_bcell.obs['cell.labels']
              .map(lambda x: {'HSC_MPP': 'HSC', 'pre-B cell': 'Pre B-cell', 'pro-B cell': 'Pro B-cell', 'Pre pro B cell': 'Pre Pro B-cell', 'B cell': 'B-cell'}.get(x, x))
              .astype('category')
)
ad_bcell = ad_bcell[ad_bcell.obs['percent.mito'] < 0.15]
ad_bcell.layers['counts'] = ad_bcell.X.copy() 
ad_bcell.layers['raw_counts'] = ad_bcell.X.copy() 
ad_bcell.layers['spliced'] = ad_bcell.X.copy() 
ad_bcell.layers['unspliced'] = ad_bcell.X.copy()
ad_bcell.raw = ad_bcell

4.2 Preprocessing

We’ll begin our analysis by running our cells through a standard preprocessing pipeline composed of QC, integration, normalization, HVG selection, dimension reduction, & clustering.

4.2.1 QC

We filter out cells with a sequencing depth of less than 1,000, then remove out genes expressed in less than 5 cells.

Code 
sc.pp.filter_cells(ad_bcell, min_counts=1000)
sc.pp.filter_genes(ad_bcell, min_cells=5)

4.2.2 HVG selection

Here we select the top 4000 most highly-variable genes (HVGs), which we’ll use as input to PCA, integration, etc.

Code 
sc.pp.highly_variable_genes(
    ad_bcell, 
    n_top_genes=4000, 
    flavor='seurat_v3', 
    layer='counts',
    subset=True
)

4.2.3 Integration with scVI

Using the scVI library, we perform celltype-aware integration across batches by training a variational autoencoder (VAE) for 150 epochs. This provides us with a 20-dimensional integrated latent space embedding that is (hopefully) free of batch effects.

Code 
scvi.settings.seed = 312
Code 
scvi.settings.num_threads = 12
scvi.model.SCVI.setup_anndata(
    ad_bcell, 
    layer='counts', 
    batch_key='fetal.ids', 
    labels_key='celltype'
)
int_model = scvi.model.SCVI(
    ad_bcell, 
    n_layers=2, 
    n_hidden=72, 
    n_latent=20, 
    gene_likelihood='nb', 
    dispersion='gene-label'
)
int_model.train(
    early_stopping=True,
    accelerator='cpu', 
    max_epochs=150
)
Code 
ad_bcell.obsm['X_scVI'] = int_model.get_latent_representation()

We plot the first two scVI components below.

Code 
sc.pl.embedding(
    ad_bcell, 
    basis='scVI', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    show=False
)
plt.gca().set_xlabel('scVI 1')
plt.gca().set_ylabel('scVI 2')
plt.show()
Figure 1: scVI embedding colored by celltype.

4.2.4 Normalization

We next depth- and log1p-normalize the mRNA counts matrix.

Code 
ad_bcell.X = sc.pp.normalize_total(ad_bcell, target_sum=1e4, inplace=False)['X']
sc.pp.log1p(ad_bcell)

4.2.5 PCA embedding

Using PCA we generate a 50-dimensional linear reduction of the normalized counts.

Code 
sc.pp.scale(ad_bcell)
sc.tl.pca(
    ad_bcell, 
    n_comps=50, 
    random_state=312, 
    use_highly_variable=True
)

Plotting the first two PCs shows some separation by celltype.

Code 
sc.pl.embedding(
    ad_bcell, 
    basis='pca', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 2: PCA embedding colored by celltype.

4.2.6 SNN graph estimation

Next, we identify the 50 nearest-neighbors (NNs) for each cell using the cosine distance in the scVI latent space.

Code 
sc.pp.neighbors(
    ad_bcell, 
    n_neighbors=50,
    n_pcs=None,  
    metric='cosine', 
    random_state=312, 
    use_rep='X_scVI'
)

4.2.7 Clustering

Using the Leiden algorithm we partition the SNN graph into clusters.

Code 
sc.tl.leiden(
    ad_bcell, 
    resolution=0.3, 
    random_state=312
)

The identified clusters roughly match our celltypes.

Code 
sc.pl.embedding(
    ad_bcell, 
    basis='pca', 
    color='leiden',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('PC 1')
plt.gca().set_ylabel('PC 2')
plt.show()
Figure 3: PCA embedding colored by Leiden cluster.

4.2.8 Moment-based imputation

Moving on, we create smoothed versions of the counts across each cell’s 30 NNs.

Code 
scv.pp.moments(
    ad_bcell, 
    n_pcs=None,
    use_rep='X_scVI', 
    n_neighbors=30
)

4.2.9 Undirected PAGA embedding

We use the PAGA algorithm to generate a graph abstraction of the relationships between celltypes.

Code 
sc.tl.paga(ad_bcell, groups='celltype')

Visualizing the results shows a pretty linear structure that matches what we expect biologically.

Code 
sc.pl.paga(
    ad_bcell, 
    fontoutline=True, 
    fontsize=12, 
    frameon=True, 
    random_state=312, 
    show=False
)
plt.gca().set_xlabel('PAGA 1')
plt.gca().set_ylabel('PAGA 2')
plt.show()
Figure 4: PAGA embedding colored by celltype. Lines are weighted by (undirected) connectivity strength.

4.2.10 UMAP embedding

We then create a nonlinear dimension reduction of the data using the UMAP algorithm.

Code 
sc.tl.umap(ad_bcell, random_state=312)

In UMAP space the transitions between celltypes are relatively clear and well-represented, though there’s a lack of connectivity between the progenitor celltypes and the mature B-cells which isn’t ideal.

Code 
sc.pl.embedding(
    ad_bcell, 
    basis='umap', 
    color='celltype',
    title='',
    frameon=True, 
    size=30, 
    alpha=0.5,
    show=False
)
plt.gca().set_xlabel('UMAP 1')
plt.gca().set_ylabel('UMAP 2')
plt.show()
Figure 5: UMAP embedding colored by celltype.

4.2.11 Force-directed graph embedding

Using the ForceAtlas2 algorithm we generate a 2-dimensional graph layout of our cells.

Code 
sc.tl.draw_graph(
    ad_bcell, 
    layout='fa', 
    random_state=312
)

The force-directed graph embedding suffers from the same issues as the UMAP embedding, namely lack of connectivity between celltypes.

Code 
sc.pl.draw_graph(
    ad_bcell, 
    color='celltype', 
    title='', 
    size=30, 
    alpha=0.5, 
    show=False, 
    frameon=True
)
plt.gca().set_xlabel('FA 1')
plt.gca().set_ylabel('FA 2')
plt.show()
Figure 6: Force-directed graph embedding colored by celltype.

4.2.12 Diffusion map embedding

Our last embedding will be generated using diffusion maps, which explicitly attempts to preserve transitional structures in the data.

Code 
sc.tl.diffmap(
    ad_bcell, 
    random_state=312, 
    n_comps=16
)
ad_bcell.obsm['X_diffmap_old'] = ad_bcell.obsm['X_diffmap']
ad_bcell.obsm['X_diffmap'] = ad_bcell.obsm['X_diffmap'][:, 1:]

The diffusion map embedding seems to correctly preserve transitions between B-cell subtypes. We’ll use this embedding going forwards.

Code 
sc.pl.embedding(
    ad_bcell, 
    basis='diffmap', 
    color='celltype',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 7: Diffusion map embedding colored by celltype.

4.2.13 Diffusion pseudotime estimation

We next compute a diffusion pseudotime estimate for each cell using the first diffusion component.

Code 
ad_bcell.uns['iroot'] = np.argmax(ad_bcell.obsm['X_diffmap'][:, 0])
sc.tl.dpt(ad_bcell, n_dcs=1)

Overall the diffusion pseudotime captures B-cell differentiation progression well. We’ll use this cellular ordering for our trajectory DE testing.

Code 
sc.pl.embedding(
    ad_bcell, 
    basis='diffmap', 
    color='dpt_pseudotime',
    title='', 
    frameon=True, 
    size=30, 
    alpha=0.5, 
    components='1,2', 
    show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 8: Diffusion map embedding colored by diffusion pseudotime.

4.3 Save preprocessed data

We save the processed AnnData object to disk.

Code 
ad_bcell.write_h5ad('../../data/bcell/ad_bcell.h5ad')

We also save the embeddings, counts matrix, metadata, etc. so that we can create a Seurat object.

Code 
ad_bcell.obs.reset_index(inplace=False).rename(columns={'index': 'cell'}).to_csv('../../data/bcell/conversion_files/cell_metadata.csv')
ad_bcell.var.reset_index(inplace=False).rename(columns={'index': 'gene'}).to_csv('../../data/bcell/conversion_files/gene_metadata.csv')
pd.DataFrame(ad_bcell.obsm['X_draw_graph_fa']).reset_index(inplace=False).to_csv('../../data/bcell/conversion_files/graph_embed.csv')
pd.DataFrame(ad_bcell.obsm['X_umap']).reset_index(inplace=False).to_csv('../../data/bcell/conversion_files/umap_embed.csv')
pd.DataFrame(ad_bcell.obsm['X_pca']).reset_index(inplace=False).to_csv('../../data/bcell/conversion_files/pca_embed.csv')
pd.DataFrame(ad_bcell.obsm['X_scVI']).reset_index(inplace=False).to_csv('../../data/bcell/conversion_files/scvi_embed.csv')
pd.DataFrame(ad_bcell.obsm['X_diffmap']).reset_index(inplace=False).to_csv('../../data/bcell/conversion_files/diffmap_embed.csv')
mmwrite('../../data/bcell/conversion_files/raw_counts.mtx', a=ad_bcell.layers['raw_counts'])

5 Analysis

5.1 Read data into R

We start our analysis by importing the preprocessed data from Python into R. Expand the code block for details.

Code 
cell_metadata <- readr::read_csv("../../data/bcell/conversion_files/cell_metadata.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1) %>% 
                 as.data.frame() %>% 
                 magrittr::set_rownames(.$cell)
gene_metadata <- readr::read_csv("../../data/bcell/conversion_files/gene_metadata.csv", 
                                 show_col_types = FALSE, 
                                 col_select = -1)
embed_PCA <- readr::read_csv("../../data/bcell/conversion_files/pca_embed.csv",
                             show_col_types = FALSE, 
                             col_select = -c(1:2)) %>% 
             as.data.frame() %>% 
             magrittr::set_colnames(paste0("PC", 1:ncol(.))) %>% 
             magrittr::set_rownames(cell_metadata$cell)
embed_PCA <- CreateDimReducObject(embeddings = as.matrix(embed_PCA),
                                  key = "PC_", 
                                  global = TRUE, 
                                  assay = "RNA")
embed_FA <- readr::read_csv("../../data/bcell/conversion_files/graph_embed.csv", 
                            show_col_types = FALSE, 
                            col_select = -c(1:2)) %>% 
            as.data.frame() %>% 
            magrittr::set_colnames(paste0("FA", 1:ncol(.))) %>% 
            magrittr::set_rownames(cell_metadata$cell)
embed_FA <- CreateDimReducObject(embeddings = as.matrix(embed_FA),
                                 key = "FA_", 
                                 global = TRUE, 
                                 assay = "RNA")
embed_UMAP <- readr::read_csv("../../data/bcell/conversion_files/umap_embed.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("UMAP", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_UMAP <- CreateDimReducObject(embeddings = as.matrix(embed_UMAP),
                                   key = "UMAP_", 
                                   global = TRUE, 
                                   assay = "RNA")
embed_diffmap <- readr::read_csv("../../data/bcell/conversion_files/diffmap_embed.csv", 
                                 show_col_types = FALSE, 
                                  col_select = -c(1:2)) %>% 
                 as.data.frame() %>% 
                 magrittr::set_colnames(paste0("DC", 1:ncol(.))) %>% 
                 magrittr::set_rownames(cell_metadata$cell)
embed_diffmap <- CreateDimReducObject(embeddings = as.matrix(embed_diffmap),
                                      key = "DC_", 
                                      global = TRUE, 
                                      assay = "RNA")
embed_scVI <- readr::read_csv("../../data/bcell/conversion_files/scvi_embed.csv", 
                              show_col_types = FALSE, 
                               col_select = -c(1:2)) %>% 
              as.data.frame() %>% 
              magrittr::set_colnames(paste0("scVI", 1:ncol(.))) %>% 
              magrittr::set_rownames(cell_metadata$cell)
embed_scVI <- CreateDimReducObject(embeddings = as.matrix(embed_scVI),
                                   key = "scVI_", 
                                   global = TRUE, 
                                   assay = "RNA")
raw_counts <- Matrix::t(Matrix::readMM("../../data/bcell/conversion_files/raw_counts.mtx"))
colnames(raw_counts) <- cell_metadata$cell
rownames(raw_counts) <- gene_metadata$gene
raw_counts <- CreateAssayObject(counts = raw_counts,
                                min.cells = 0,
                                min.features = 0, 
                                key = "rna")

We can now create a Seurat object, to which we add our various embeddings and metadata. After adding all the various bits and pieces, we normalize the raw counts and re-identify HVGs.

Code 
seu_bcell <- CreateSeuratObject(raw_counts, 
                                assay = "RNA",
                                meta.data = cell_metadata, 
                                project = "bcell", 
                                min.cells = 0, 
                                min.features = 0)
seu_bcell@meta.data <- mutate(seu_bcell@meta.data, 
                              celltype = factor(celltype, levels = c("HSC", 
                                                                     "Pre Pro B-cell", 
                                                                     "Pro B-cell", 
                                                                     "Pre B-cell", 
                                                                     "B-cell")))
seu_bcell@reductions$pca <- embed_PCA
seu_bcell@reductions$diffmap <- embed_diffmap
seu_bcell@reductions$fa <- embed_FA
seu_bcell@reductions$umap <- embed_UMAP
seu_bcell@reductions$scvi <- embed_scVI
seu_bcell <- NormalizeData(seu_bcell, verbose = FALSE) %>% 
             FindVariableFeatures(selection.method = "vst", verbose = FALSE)

5.2 scLANE DE analysis

We use the scLANE package to perform trajectory DE testing across pseudotime for a set of 2,000 candidate genes. Since our dataset is of moderate size and composed of samples from multiple subjects we’ll utilize the GEE mode of scLANE.

Code 
candidate_genes <- chooseCandidateGenes(seu_bcell,
                                        id.vec = seu_bcell$fetal.ids, 
                                        n.desired.genes = 2000L)
seu_bcell_sort <- sortObservations(seu_bcell, 
                                   pt = seu_bcell$dpt_pseudotime, 
                                   id.vec = seu_bcell$fetal.ids)
Warning in sortObservations(seu_bcell, pt = seu_bcell$dpt_pseudotime, id.vec =
seu_bcell$fetal.ids): Ordering a Seurat object requires conversion to
SingleCellExperiment, and some information might be lost.
Code 
cell_offset <- createCellOffset(seu_bcell_sort)
pt_df <- data.frame(DPT = seu_bcell_sort$dpt_pseudotime)
scLANE_models <- testDynamic(seu_bcell_sort, 
                             pt = pt_df, 
                             genes = candidate_genes, 
                             size.factor.offset = cell_offset, 
                             is.gee = TRUE, 
                             id.vec = seu_bcell_sort$fetal.ids, 
                             n.cores = 24L,
                             verbose = FALSE)
scLANE testing in GEE mode completed for 2000 genes across 1 lineage in 2.125 hours

5.3 tradeSeq DE analysis

Continuing on, we perform subject-naive TDE testing via tradeSeq.

Code 
ts_start <- Sys.time()
bioc_par <- BiocParallel::MulticoreParam(workers = 24L, RNGseed = 312)
RNA_counts <- as.matrix(seu_bcell@assays$RNA$counts)[candidate_genes, ]
cell_offset_2 <- createCellOffset(seu_bcell)
pt_df_2 <- data.frame(DPT = seu_bcell$dpt_pseudotime)
k_eval <- evaluateK(RNA_counts, 
                    pseudotime = pt_df_2, 
                    cellWeights = matrix(rep(1, nrow(pt_df_2)), ncol = 1), 
                    offset = log(1 / cell_offset_2), 
                    k = 3:10, 
                    plot = FALSE, 
                    nGenes = 500, 
                    verbose = FALSE, 
                    parallel = TRUE, 
                    BPPARAM = bioc_par)
best_k <- c(3:10)[which.min(abs(colMeans(k_eval - rowMeans(k_eval))))]  # choose k w/ lowest MAD from mean AIC 
ts_models <- fitGAM(RNA_counts, 
                    pseudotime = pt_df_2, 
                    cellWeights = matrix(rep(1, nrow(pt_df_2)), ncol = 1), 
                    offset = log(1 / cell_offset_2), 
                    nknots = best_k, 
                    sce = FALSE, 
                    parallel = TRUE, 
                    verbose = FALSE, 
                    BPPARAM = bioc_par)
BiocParallel::bpstop(bioc_par)
ts_end <- Sys.time()
ts_diff <- ts_end - ts_start

The runtime for tradeSeq with 5 knots per gene is:

Code 
ts_diff
Time difference of 18.33945 mins

5.4 Lamian DE analysis

Next, we perform subject-aware trajectory DE testing with Lamian using their default settings.

Code 
lamian_start <- Sys.time()
cell_anno <- select(seu_bcell@meta.data, 
                    Cell = cell, 
                    Sample = fetal.ids) %>% 
             mutate(across(everything(), as.character)) %>% 
             magrittr::set_rownames(NULL)
cell_pt <- as.integer(rank(seu_bcell$dpt_pseudotime))
names(cell_pt) <- cell_anno$Cell
samp_design <- data.frame(intercept = rep(1, length(unique(seu_bcell$fetal.ids)))) %>%
               magrittr::set_rownames(unique(seu_bcell$fetal.ids)) %>%
               as.matrix()
lamian_models <- lamian_test(expr = as.matrix(seu_bcell@assays$RNA$data)[candidate_genes, ],
                             cellanno = cell_anno,
                             pseudotime = cell_pt,
                             design = samp_design,
                             test.type = "time", 
                             test.method = "permutation", 
                             permuiter = 100, 
                             verbose.output = FALSE, 
                             ncores = 24L)
lamian_end <- Sys.time()
lamian_diff <- lamian_end - lamian_start

The runtime for Lamian is:

Code 
lamian_diff
Time difference of 18.03005 mins

5.5 Cell fate analysis

Lastly, we’ll perform an analysis of cell fate specification using the CellRank Python package.

5.5.1 CytoTRACE kernel

We begin by computing a CytoTRACE kernel; CytoTRACE uses the number of genes expressed in a cell as a proxy measurement for differentiation potential, with the assumption being that cells expressing more genes are more likely to be progenitor celltypes & vice versa.

Code 
ctk = CytoTRACEKernel(ad_bcell).compute_cytotrace()

We plot the CytoTRACE scores & pseudotime below.

Code 
sc.pl.embedding(
    ad_bcell,
    basis='diffmap',
    color=['ct_score', 'ct_pseudotime'],
    color_map='magma', 
    show=False, 
    size=30, 
    alpha=0.5, 
    title=['CytoTRACE score', 'CytoTRACE Pseudotime']
)
plt.show()
Figure 9: Diffusion map embedding colored by CytoTRACE score and pseudotime. Higher scores indicate higher differentiation potential.

Next, we estimate a cell-cell transition probability matrix.

Code 
ctk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the projection of that matrix on our diffusion map embedding shows messy but overall correct directionality.

Code 
ctk.plot_projection(
  basis='diffmap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.5, 
  linewidth=2, 
  components='0,1',
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 10: Diffusion map embedding streamline plot showing differentiation directions as inferred from CytoTRACE scores.

5.5.2 Pseudotime kernel

Next, we use the diffusion pseudotime estimates to create a pseudotime-based kernel and estimate another cell-cell transition probability matrix.

Code 
pk = PseudotimeKernel(ad_bcell, time_key='dpt_pseudotime')
pk.compute_transition_matrix(threshold_scheme='soft', show_progress_bar=False)

Plotting the projection of the pseudotime kernel’s transition probability matrix shows a more reasonable portrayal of differentiation directionality.

Code 
pk.plot_projection(
  basis='diffmap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  components='0,1',
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 11: Diffusion map embedding overlaid with streamline arrows derived from Slingshot pseudotime.

5.5.3 Combined kernel

We combine the two kernel using unequal weighting, with the pseudotime kernel being given priority due to its better overall capture of differentiation.

Code 
ck = 0.3 * ctk + 0.7 * pk

We visualize the combined kernel’s projection below. It seems to reliably portray the direction of differentiation in our dataset.

Code 
ck.plot_projection(
  basis='diffmap', 
  recompute=True, 
  color='celltype',
  title='',
  legend_loc='right margin', 
  size=30, 
  alpha=0.75, 
  linewidth=2, 
  components='0,1',
  frameon=True, 
  show=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 12: Diffusion map embedding overlaid with streamline arrows derived from the combined kernel.

We can also simulate random walks on the high-dimensional cell-cell graph starting from the HSC population. We see that nearly all of the random walks terminate in the mature B-cell population.

Code 
ck.plot_random_walks(
    n_sims=200,
    start_ixs={'celltype': 'HSC'},
    basis='diffmap',
    color='celltype',
    legend_loc='right margin',
    seed=312, 
    frameon=True, 
    title='', 
    linewidth=0.5, 
    linealpha=0.25, 
    size=30, 
    alpha=0.5, 
    components='0,1',
    ixs_legend_loc='upper', 
    show_progress_bar=False
)
plt.gca().set_xlabel('DC 1')
plt.gca().set_ylabel('DC 2')
plt.show()
Figure 13: Diffusion map embedding overlaid with random walks along the cell-cell graph derived from the combined kernel. Starting points are highlighted in black and ending points in yellow.

6 Save data

Finally, we save our Seurat object, the output from Slingshot, and the models from scLANE and tradeSeq.

Code 
saveRDS(seu_bcell, file = "../../data/bcell/seu_bcell.Rds")
saveRDS(scLANE_models, file = "../../data/bcell/scLANE_models.Rds")
saveRDS(ts_models, file = "../../data/bcell/ts_models.Rds")
saveRDS(lamian_models, file = "../../data/bcell/lamian_models.Rds")

7 Session info

Code 
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.3 (2024-02-29)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/New_York
 date     2024-10-06
 pandoc   2.9.2.1 @ /usr/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package              * version    date (UTC) lib source
 abind                  1.4-8      2024-09-12 [2] CRAN (R 4.3.3)
 AnnotationDbi          1.62.2     2023-07-02 [2] Bioconductor
 backports              1.5.0      2024-05-23 [1] CRAN (R 4.2.3)
 bigassertr             0.1.6      2023-01-10 [1] CRAN (R 4.3.2)
 bigparallelr           0.3.2      2021-10-02 [1] CRAN (R 4.3.2)
 bigstatsr              1.5.12     2022-10-14 [1] CRAN (R 4.3.2)
 Biobase              * 2.60.0     2023-04-25 [2] Bioconductor
 BiocFileCache          2.13.1     2024-10-05 [1] Github (Bioconductor/BiocFileCache@30abbe4)
 BiocGenerics         * 0.46.0     2023-04-25 [2] Bioconductor
 BiocParallel           1.34.2     2023-05-22 [2] Bioconductor
 biomaRt              * 2.56.1     2023-06-09 [1] Bioconductor
 Biostrings             2.68.1     2023-05-16 [2] Bioconductor
 bit                    4.5.0      2024-09-20 [2] CRAN (R 4.3.3)
 bit64                  4.5.2      2024-09-22 [2] CRAN (R 4.3.3)
 bitops                 1.0-8      2024-07-29 [2] CRAN (R 4.3.3)
 blob                   1.2.4      2023-03-17 [2] CRAN (R 4.3.1)
 boot                   1.3-31     2024-08-28 [4] CRAN (R 4.3.3)
 broom                  1.0.6      2024-05-17 [1] CRAN (R 4.2.3)
 broom.mixed            0.2.9.5    2024-04-01 [1] CRAN (R 4.2.3)
 cachem                 1.1.0      2024-05-16 [2] CRAN (R 4.3.3)
 caTools                1.18.3     2024-09-04 [2] CRAN (R 4.3.3)
 circlize               0.4.16     2024-02-20 [1] CRAN (R 4.3.2)
 cli                    3.6.2      2023-12-11 [1] CRAN (R 4.2.3)
 clue                   0.3-65     2023-09-23 [1] CRAN (R 4.3.2)
 cluster                2.1.6      2023-12-01 [4] CRAN (R 4.3.2)
 coda                   0.19-4.1   2024-01-31 [2] CRAN (R 4.3.2)
 codetools              0.2-20     2024-03-31 [4] CRAN (R 4.3.3)
 colorspace             2.1-1      2024-07-26 [2] CRAN (R 4.3.3)
 combinat               0.0-8      2012-10-29 [2] CRAN (R 4.3.1)
 ComplexHeatmap         2.16.0     2023-04-25 [1] Bioconductor
 cowplot                1.1.3      2024-01-22 [2] CRAN (R 4.3.3)
 crayon                 1.5.3      2024-06-20 [2] CRAN (R 4.3.3)
 curl                   5.2.3      2024-09-20 [2] CRAN (R 4.3.3)
 data.table             1.16.0     2024-08-27 [2] CRAN (R 4.3.3)
 DBI                    1.2.3      2024-06-02 [2] CRAN (R 4.3.3)
 dbplyr                 2.5.0      2024-03-19 [1] CRAN (R 4.3.3)
 DelayedArray           0.26.7     2023-07-28 [2] Bioconductor
 deldir                 2.0-4      2024-02-28 [2] CRAN (R 4.3.3)
 devtools               2.4.5      2022-10-11 [1] CRAN (R 4.2.2)
 digest                 0.6.35     2024-03-11 [1] CRAN (R 4.2.3)
 doParallel             1.0.17     2022-02-07 [1] CRAN (R 4.2.3)
 doSNOW                 1.0.20     2022-02-04 [1] CRAN (R 4.3.2)
 dotCall64              1.1-1      2023-11-28 [2] CRAN (R 4.3.1)
 dplyr                * 1.1.4      2023-11-17 [1] CRAN (R 4.2.3)
 edgeR                  3.42.4     2023-05-31 [2] Bioconductor
 ellipsis               0.3.2      2021-04-29 [2] CRAN (R 4.3.1)
 emmeans                1.10.2     2024-05-20 [1] CRAN (R 4.3.3)
 estimability           1.5.1      2024-05-12 [1] CRAN (R 4.3.3)
 evaluate               0.24.0     2024-06-10 [1] CRAN (R 4.3.3)
 fansi                  1.0.6      2023-12-08 [2] CRAN (R 4.3.1)
 farver                 2.1.2      2024-05-13 [1] CRAN (R 4.2.3)
 fastDummies            1.7.4      2024-08-16 [2] CRAN (R 4.3.3)
 fastICA                1.2-4      2023-11-27 [1] CRAN (R 4.3.2)
 fastmap                1.2.0      2024-05-15 [2] CRAN (R 4.3.3)
 ff                     4.0.12     2024-01-12 [1] CRAN (R 4.3.2)
 filelock               1.0.3      2023-12-11 [2] CRAN (R 4.3.3)
 fitdistrplus           1.2-1      2024-07-12 [2] CRAN (R 4.3.3)
 flock                  0.7        2016-11-12 [1] CRAN (R 4.3.2)
 forcats                1.0.0      2023-01-29 [1] CRAN (R 4.3.2)
 foreach                1.5.2      2022-02-02 [1] CRAN (R 4.2.3)
 fs                     1.6.4      2024-04-25 [2] CRAN (R 4.3.3)
 furrr                  0.3.1      2022-08-15 [1] CRAN (R 4.3.2)
 future                 1.33.2     2024-03-26 [1] CRAN (R 4.2.3)
 future.apply           1.11.2     2024-03-28 [2] CRAN (R 4.3.3)
 gamlss                 5.4-22     2024-03-20 [1] CRAN (R 4.3.2)
 gamlss.data            6.0-6      2024-03-14 [1] CRAN (R 4.3.2)
 gamlss.dist            6.1-1      2023-08-23 [1] CRAN (R 4.3.2)
 geeM                   0.10.1     2018-06-18 [1] CRAN (R 4.3.2)
 generics               0.1.3      2022-07-05 [1] CRAN (R 4.2.0)
 GenomeInfoDb         * 1.36.4     2023-10-02 [2] Bioconductor
 GenomeInfoDbData       1.2.10     2023-07-26 [2] Bioconductor
 GenomicRanges        * 1.52.1     2023-10-08 [2] Bioconductor
 GetoptLong             1.0.5      2020-12-15 [1] CRAN (R 4.3.2)
 ggplot2              * 3.5.1      2024-04-23 [1] CRAN (R 4.2.3)
 ggrepel                0.9.2      2022-11-06 [1] CRAN (R 4.2.2)
 ggridges               0.5.6      2024-01-23 [2] CRAN (R 4.3.3)
 glm2                 * 1.2.1      2018-08-11 [1] CRAN (R 4.3.2)
 glmmTMB                1.1.9      2024-03-20 [1] CRAN (R 4.3.3)
 GlobalOptions          0.1.2      2020-06-10 [1] CRAN (R 4.3.2)
 globals                0.16.3     2024-03-08 [2] CRAN (R 4.3.3)
 glue                   1.8.0      2024-09-30 [2] CRAN (R 4.3.3)
 GO.db                  3.17.0     2023-08-11 [2] Bioconductor
 goftest                1.2-3      2021-10-07 [2] CRAN (R 4.3.1)
 gplots                 3.1.3.1    2024-02-02 [2] CRAN (R 4.3.3)
 graph                  1.78.0     2023-04-25 [1] Bioconductor
 gridExtra              2.3        2017-09-09 [2] CRAN (R 4.3.1)
 gtable                 0.3.5      2024-04-22 [2] CRAN (R 4.3.3)
 gtools                 3.9.5      2023-11-20 [2] CRAN (R 4.3.1)
 hms                    1.1.3      2023-03-21 [2] CRAN (R 4.3.1)
 htmltools              0.5.8.1    2024-04-04 [1] CRAN (R 4.2.3)
 htmlwidgets            1.6.4      2023-12-06 [2] CRAN (R 4.3.3)
 httpuv                 1.6.15     2024-03-26 [2] CRAN (R 4.3.3)
 httr                   1.4.7      2023-08-15 [2] CRAN (R 4.3.3)
 ica                    1.0-3      2022-07-08 [2] CRAN (R 4.3.1)
 igraph                 2.0.3      2024-03-13 [1] CRAN (R 4.2.3)
 IRanges              * 2.34.1     2023-06-22 [2] Bioconductor
 irlba                  2.3.5.1    2022-10-03 [2] CRAN (R 4.3.1)
 iterators              1.0.14     2022-02-05 [2] CRAN (R 4.3.1)
 jsonlite               1.8.8      2023-12-04 [1] CRAN (R 4.2.3)
 KEGGREST               1.40.1     2023-09-29 [2] Bioconductor
 KernSmooth             2.23-24    2024-05-17 [4] CRAN (R 4.3.3)
 knitr                  1.48       2024-07-07 [2] CRAN (R 4.3.3)
 Lamian               * 0.99.2     2024-03-20 [1] Github (Winnie09/Lamian@7d5a8ff)
 later                  1.3.2      2023-12-06 [2] CRAN (R 4.3.3)
 lattice                0.22-6     2024-03-20 [4] CRAN (R 4.3.3)
 lazyeval               0.2.2      2019-03-15 [2] CRAN (R 4.3.1)
 leiden                 0.4.3.1    2023-11-17 [2] CRAN (R 4.3.1)
 lifecycle              1.0.4      2023-11-07 [2] CRAN (R 4.3.1)
 limma                  3.56.2     2023-06-04 [2] Bioconductor
 listenv                0.9.1      2024-01-29 [2] CRAN (R 4.3.3)
 lme4                   1.1-35.5   2024-07-03 [2] CRAN (R 4.3.3)
 lmtest                 0.9-40     2022-03-21 [2] CRAN (R 4.3.1)
 locfit                 1.5-9.10   2024-06-24 [2] CRAN (R 4.3.3)
 magrittr             * 2.0.3      2022-03-30 [1] CRAN (R 4.2.0)
 MASS                   7.3-60     2023-05-04 [4] CRAN (R 4.3.1)
 Matrix                 1.6-5      2024-01-11 [1] CRAN (R 4.2.3)
 matrixcalc             1.0-6      2022-09-14 [1] CRAN (R 4.3.2)
 MatrixGenerics       * 1.12.3     2023-07-30 [2] Bioconductor
 matrixStats          * 1.4.1      2024-09-08 [2] CRAN (R 4.3.3)
 mclust                 6.1.1      2024-04-29 [2] CRAN (R 4.3.3)
 memoise                2.0.1      2021-11-26 [2] CRAN (R 4.3.1)
 mgcv                   1.9-1      2023-12-21 [4] CRAN (R 4.3.2)
 mime                   0.12       2021-09-28 [2] CRAN (R 4.3.1)
 miniUI                 0.1.1.1    2018-05-18 [2] CRAN (R 4.3.1)
 minqa                  1.2.7      2024-05-20 [1] CRAN (R 4.2.3)
 munsell                0.5.1      2024-04-01 [2] CRAN (R 4.3.3)
 mvtnorm                1.3-1      2024-09-03 [2] CRAN (R 4.3.3)
 nlme                   3.1-166    2024-08-14 [4] CRAN (R 4.3.3)
 nloptr                 2.1.1      2024-06-25 [2] CRAN (R 4.3.3)
 numDeriv               2016.8-1.1 2019-06-06 [2] CRAN (R 4.3.1)
 org.Hs.eg.db           3.17.0     2024-01-31 [2] Bioconductor
 paletteer              1.6.0      2024-01-21 [1] CRAN (R 4.2.3)
 parallelly             1.37.1     2024-02-29 [1] CRAN (R 4.2.3)
 patchwork            * 1.2.0      2024-01-08 [1] CRAN (R 4.2.3)
 pbapply                1.7-2      2023-06-27 [2] CRAN (R 4.3.1)
 pheatmap               1.0.12     2019-01-04 [1] CRAN (R 4.3.3)
 pillar                 1.9.0      2023-03-22 [2] CRAN (R 4.3.1)
 pkgbuild               1.4.4      2024-03-17 [2] CRAN (R 4.3.3)
 pkgconfig              2.0.3      2019-09-22 [2] CRAN (R 4.3.1)
 pkgload                1.4.0      2024-06-28 [2] CRAN (R 4.3.3)
 plotly                 4.10.4     2024-01-13 [2] CRAN (R 4.3.3)
 plyr                   1.8.9      2023-10-02 [2] CRAN (R 4.3.3)
 png                    0.1-8      2022-11-29 [2] CRAN (R 4.3.1)
 polyclip               1.10-7     2024-07-23 [2] CRAN (R 4.3.3)
 prettyunits            1.2.0      2023-09-24 [2] CRAN (R 4.3.3)
 princurve            * 2.1.6      2021-01-18 [1] CRAN (R 4.3.2)
 prismatic              1.1.2      2024-04-10 [1] CRAN (R 4.2.3)
 profvis                0.4.0      2024-09-20 [2] CRAN (R 4.3.3)
 progress               1.2.3      2023-12-06 [2] CRAN (R 4.3.3)
 progressr              0.14.0     2023-08-10 [2] CRAN (R 4.3.2)
 promises               1.3.0      2024-04-05 [2] CRAN (R 4.3.3)
 ps                     1.8.0      2024-09-12 [2] CRAN (R 4.3.3)
 purrr                  1.0.2      2023-08-10 [1] CRAN (R 4.2.3)
 R6                     2.5.1      2021-08-19 [2] CRAN (R 4.3.1)
 RANN                   2.6.2      2024-08-25 [2] CRAN (R 4.3.3)
 rappdirs               0.3.3      2021-01-31 [2] CRAN (R 4.3.1)
 RColorBrewer           1.1-3      2022-04-03 [2] CRAN (R 4.3.1)
 Rcpp                   1.0.13     2024-07-17 [2] CRAN (R 4.3.3)
 RcppAnnoy              0.0.22     2024-01-23 [2] CRAN (R 4.3.3)
 RcppEigen              0.3.4.0.0  2024-02-28 [1] CRAN (R 4.2.3)
 RcppHNSW               0.6.0      2024-02-04 [2] CRAN (R 4.3.3)
 RCurl                  1.98-1.16  2024-07-11 [2] CRAN (R 4.3.3)
 readr                  2.1.5      2024-01-10 [2] CRAN (R 4.3.3)
 rematch2               2.1.2      2020-05-01 [2] CRAN (R 4.3.1)
 remotes                2.5.0      2024-03-17 [2] CRAN (R 4.3.3)
 reshape2               1.4.4      2020-04-09 [2] CRAN (R 4.3.1)
 reticulate           * 1.35.0     2024-01-31 [1] CRAN (R 4.2.3)
 rhdf5                  2.44.0     2023-04-25 [2] Bioconductor
 rhdf5filters           1.12.1     2023-04-30 [2] Bioconductor
 Rhdf5lib               1.22.1     2023-09-10 [2] Bioconductor
 rjson                  0.2.23     2024-09-16 [2] CRAN (R 4.3.3)
 rlang                  1.1.4      2024-06-04 [1] CRAN (R 4.3.3)
 rmarkdown              2.27       2024-05-17 [1] CRAN (R 4.3.3)
 rmio                   0.4.0      2022-02-17 [1] CRAN (R 4.3.2)
 ROCR                   1.0-11     2020-05-02 [2] CRAN (R 4.3.1)
 RSpectra               0.16-2     2024-07-18 [2] CRAN (R 4.3.3)
 RSQLite                2.3.7      2024-05-27 [2] CRAN (R 4.3.3)
 Rtsne                  0.17       2023-12-07 [2] CRAN (R 4.3.3)
 S4Arrays               1.2.1      2024-03-04 [1] Bioconductor 3.18 (R 4.3.2)
 S4Vectors            * 0.38.2     2023-09-22 [2] Bioconductor
 scales                 1.3.0      2023-11-28 [1] CRAN (R 4.2.3)
 scattermore            1.2        2023-06-12 [2] CRAN (R 4.3.1)
 scLANE               * 0.8.4      2024-10-05 [1] Github (jr-leary7/scLANE@d1eb84d)
 sctransform            0.4.1      2023-10-19 [1] CRAN (R 4.2.3)
 sessioninfo            1.2.2      2021-12-06 [2] CRAN (R 4.3.1)
 Seurat               * 5.1.0      2024-05-10 [1] CRAN (R 4.3.3)
 SeuratObject         * 5.0.2      2024-05-08 [1] CRAN (R 4.3.3)
 shape                  1.4.6.1    2024-02-23 [2] CRAN (R 4.3.3)
 shiny                  1.9.1      2024-08-01 [2] CRAN (R 4.3.3)
 SingleCellExperiment * 1.22.0     2023-04-25 [2] Bioconductor
 slingshot            * 2.6.0      2022-11-01 [1] Bioconductor
 snow                   0.4-4      2021-10-27 [2] CRAN (R 4.3.1)
 sp                   * 2.1-4      2024-04-30 [3] CRAN (R 4.3.3)
 spam                   2.10-0     2023-10-23 [2] CRAN (R 4.3.1)
 SparseM                1.84-2     2024-07-17 [2] CRAN (R 4.3.3)
 spatstat.data          3.1-2      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.explore       3.3-2      2024-08-21 [2] CRAN (R 4.3.3)
 spatstat.geom          3.3-3      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.random        3.3-2      2024-09-18 [2] CRAN (R 4.3.3)
 spatstat.sparse        3.1-0      2024-06-21 [2] CRAN (R 4.3.3)
 spatstat.univar        3.0-1      2024-09-05 [2] CRAN (R 4.3.3)
 spatstat.utils         3.1-0      2024-08-17 [2] CRAN (R 4.3.3)
 stringi                1.8.4      2024-05-06 [2] CRAN (R 4.3.3)
 stringr                1.5.1      2023-11-14 [2] CRAN (R 4.3.2)
 SummarizedExperiment * 1.30.2     2023-06-06 [2] Bioconductor
 survival               3.7-0      2024-06-05 [4] CRAN (R 4.3.3)
 tensor                 1.5        2012-05-05 [2] CRAN (R 4.3.1)
 tibble                 3.2.1      2023-03-20 [2] CRAN (R 4.3.1)
 tidyr                  1.3.1      2024-01-24 [1] CRAN (R 4.2.3)
 tidyselect             1.2.1      2024-03-11 [2] CRAN (R 4.3.3)
 TMB                    1.9.11     2024-04-03 [1] CRAN (R 4.3.3)
 topGO                  2.52.0     2023-04-25 [1] Bioconductor
 tradeSeq             * 1.14.0     2023-04-25 [1] Bioconductor
 TrajectoryUtils      * 1.6.0      2022-11-01 [1] Bioconductor
 TSCAN                  1.36.0     2022-11-01 [1] Bioconductor
 tzdb                   0.4.0      2023-05-12 [2] CRAN (R 4.3.1)
 urlchecker             1.0.1      2021-11-30 [2] CRAN (R 4.3.1)
 usethis                2.2.3      2024-02-19 [1] CRAN (R 4.2.3)
 utf8                   1.2.4      2023-10-22 [2] CRAN (R 4.3.1)
 uwot                   0.2.2      2024-04-21 [2] CRAN (R 4.3.3)
 vctrs                  0.6.5      2023-12-01 [2] CRAN (R 4.3.1)
 viridis                0.6.5      2024-01-29 [2] CRAN (R 4.3.3)
 viridisLite            0.4.2      2023-05-02 [2] CRAN (R 4.3.1)
 vroom                  1.6.5      2023-12-05 [2] CRAN (R 4.3.3)
 withr                  3.0.0      2024-01-16 [1] CRAN (R 4.2.3)
 xfun                   0.44       2024-05-15 [1] CRAN (R 4.3.3)
 XML                    3.99-0.17  2024-06-25 [2] CRAN (R 4.3.3)
 xml2                   1.3.6      2023-12-04 [2] CRAN (R 4.3.3)
 xtable                 1.8-4      2019-04-21 [2] CRAN (R 4.3.1)
 XVector                0.40.0     2023-04-25 [2] Bioconductor
 yaml                   2.3.10     2024-07-26 [2] CRAN (R 4.3.3)
 zlibbioc               1.46.0     2023-04-25 [2] Bioconductor
 zoo                    1.8-12     2023-04-13 [2] CRAN (R 4.3.1)

 [1] /home/j.leary/r_packages_default
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /blue/rbacher/j.leary/py_envs/scLANE_env2/bin/python
 libpython:      /blue/rbacher/j.leary/py_envs/scLANE_env2/lib/libpython3.10.so
 pythonhome:     /blue/rbacher/j.leary/py_envs/scLANE_env2:/blue/rbacher/j.leary/py_envs/scLANE_env2
 version:        3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
 numpy:          /home/j.leary/.local/lib/python3.10/site-packages/numpy
 numpy_version:  1.24.2
 
 NOTE: Python version was forced by use_python() function

──────────────────────────────────────────────────────────────────────────────